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Abstract 

An analytical solution is developed for three-dimensional flow towards a partially penetrating 
large-diameter well in an unconfined aquifer bounded below by a leaky aquitard of finite or 
semi- infinite extent. The analytical solution is derived using Laplace and Hankel transforms, 
then inverted numerically. Existing solutions for flow in leaky unconfined aquifers neglect 
the unsaturated zone following an assumption of instantaneous drainage due to Neuman. 
We extend the theory of leakage in unconfined aquifers by (1) including water flow and 
storage in the unsaturated zone above the water table, and (2) allowing the finite-diameter 
pumping well to partially penetrate the aquifer. The investigation of model-predicted results 
shows that aquitard leakage leads to significant departure from the unconfined solution 
without leakage. The investigation of dimensionless time-drawdown relationships shows that 
the aquitard drawdown also depends on unsaturated zone properties and the pumping-well 
wellbore storage effects. 

Keywords: Unconfined aquifer, Aquitard, Leakage, Wellbore storage, unsaturated zone, 
Laplace-Hankel transform, Dealyed piezometer response 

1. Introduction 

The assumption that the water flow and storage in the unsaturated zone is insignificant for 
unconfined aquifer tests was first questioned by Nwankwor et al. [IJ and later by Akindunni 
and Gillham [2] based upon analysis of data collected during pumping tests in Borden, 
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5 Ontario Canada. Analyzing the collected tensiometer data and soil moisture measurements, 

6 the authors concluded that the proper inclusion of unsaturated zone in analytical models 

7 used for pumping test analysis would lead to improved estimates of aquifer specific yield. 

8 Several analytical solutions were developed that account for the unsaturated zone flow to a 

9 pumping well in an unconfined aquifer, taking into account the unsaturated zone HI [S] . 

10 These models consider the unsaturated zone effects by coupling the governing flow equations 

11 at the water table; the saturated zone governed by the diffusion equation and the vadose 

12 zone governed by the linearized unsaturated zone Richards equation, using the linearization 

13 of Kroszynski and Dagan [6] . These models considered the limiting case where the pumping 

14 well has zero radius. For detailed discussion regarding the fundamental differences between 

15 these three models readers are directed to Mishra and Neuman [3]. 

16 Drawdown due to pumping a large-diameter (e.g., water supply) well in an unconfined 

17 aquifer is affected by wellbore storage [7]. Narasimhan and Zhu [8] used a numerical model 

18 to demonstrate that early time drawdown in an unconfined aquifer tends to be dominated 

19 by wellbore storage effects. Mishra and Neuman [U] developed an analytical unconfined so- 

20 lution, which considers both pumping-well wellbore storage capacity, and three-dimensional 

21 axi-symmetrical unsaturated zone flow. They represented unsaturated zone constitutive 

22 properties using exponential models, which result in governing equations that are mathe- 

23 matically tractable, while being sufficiently flexible to be fit to other widely used constitutive 

24 models [101 [HI [iSl [13] . However, Mishra and Neuman [9] considered the unconfined aquifer 

25 to be resting on an impermeable boundary and therefore did not account for the potential 

26 effects of leakage from an underlying formation (e.g., an aquitard or fractured bedrock). 

27 The classical theory of leakage for confined aquifers was originally developed by Hantush 

28 [T3] assuming steady-state vertical flow in overlying and underlying aquitards and horizon- 

29 tal flow in the pumped aquifer. Hantush [T3] later modified the theory of confined leaky 

30 aquifers to include transient vertical aquitard flow, giving asymptotic expressions for early 

31 and late times. Neuman and Witherspoon [161 HZj developed a more complete analytical 

32 solution for the more general multiple aquifer flow problem, but did not consider general 
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33 three-dimensional aquitard flow. 

34 Yatov [18] first investigated the effect of leakage from underlying strata on unconfined 

35 aquifer flow. He use the model of Boulton [19j to account for the water table and considered 

36 only vertical flow in aquitard. Ehlig and Halepaska [20] investigated leaky-unconfined flow 

37 through a finite-difference simulation, which coupled the Boulton [H] and Hantush 

38 models to account for leakage across the aquifer-aquitard boundary. Zlotnik and Zhan [2T] 

39 developed an analytical solution for the flow towards a fully penetrating zero-radius well 

40 in a coupled unconfined aquiferaquitard system where both the unsaturated zone and the 

41 horizontal aquitard flow are neglected. Both Zhan and Bian [22j and Zlotnik and Zhan pTj 

42 developed analytical and semi-analytical solutions for leakage due to pumping, building on 

43 the works of Hantush [TJ] and Butler Jr and Tsou [23j . Both Zhan and Bian [22] and Zlotnik 

44 and Zhan [21] neglect horizontal flow in aquitards. Purely vertical aquitard flow was justified 

45 for limiting aquifer /aquitard hydraulic conductivity contrasts by Neuman and Witherspoon 

46 [TT]. Both Zhan and Bian [22] and Zlotnik and Zhan [21] only consider a vertically unbounded 

47 aquitard. Malama et al. [21] developed a solution for three-dimensional aquitard flow in a 

48 finite thickness aquitard, but considered the zero-radius pumping well to be fully penetrating 

49 and ignored the flow in unsaturated zone. Here, we develop a more general leaky-unconfined 

50 aquifer solution by considering a partially penetrating large-diameter well and including the 

51 effects of unsaturated zone flow following Mishra and Neuman |9|. The solution is used 

52 to investigate the effect of an aquitard on drawdown in overlying unconfined aquifer. We 

53 conclude by investigating the effects of wellbore storage capacity and the unsaturated zone 

54 on drawdown observed in the aquitard. 

55 2. Leaky-Unconfined Theory 

56 2.1. Statement of Problem 

57 We consider an infinite radial compressible unconfined aquifer above a finitely thick 

58 aquitard (Figure [T]). The aquifer and aquitard are each spatially uniform, homogeneous and 

59 anisotropic, with constant specific storage Sg and Sgi, respectively (a subscript 1 indicates 
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60 aquitard related properties). The aquifer has a fixed ratio Kj:, = K^/ Kr between vertical and 

61 horizontal saturated hydraulic conductivities, Kz and Kr, respectively. The aquitard vertical 

62 and horizontal hydraulic conductivities are Kzi and K^i- The aquifer is fully saturated 

63 below an initially horizontal water table at elevation z = b. The water table is defined as a 

64 ip = isobar where ip is pressure head. A saturated capillary fringe at non-positive pressure 

65 V'a < < extends from the water table to the ip = i^a isobar; 'j/'a < is the pressure 

66 head required for air to enter a saturated medium. The saturated hydraulic system (aquifer 

67 and aquitard) is at uniform initial hydraulic head = b + ipa before pumping. At time 

68 t = 0, pumping begins at a constant volumetric fiowrate Q from a well of finite radius 

69 and wellbore storage coefficient (volume of water released from storage in the pumping 

70 well per unit drawdown in the well casing). The pumping well is completed across the 

71 aquifer between depths / and d below the aquifer top. Under these conditions the drawdown 

72 s (r, z,t) = h (r, z,0) — h (r, z, t) in the saturated zone is governed by the diffusion equation 

1 d / ds\ d'^s ^ ds 

^r-JT + ^-TT^ = ^^1^ r>r^ 0<z<b, 1 

r or \ or J Oz^ Ot 

73 along with far-field boundary condition 

s(cx),2,t) = 0, (2) 

74 the no-fiow condition at the portion of the well casing that is not open to the aquifer 

ds\ 

r-T^ =0 0<z<b-lb-d<z<b, (3) 

or I 

75 and the wellbore storage mass-balance expression 

2T^Kr {I - d) ( ) - ( ^ ) =-Q b - I < z < b - d. (4) 



dr J _ \dt 

76 Flux is assumed constant across the well screen (see Zhan and Zlotnik [25j for a discussion 

77 of this assumption's validity). The corresponding linearized unsaturated fiow equations [S] 

78 are 



, , ,1 d f da\ d f . , , d(T\ ^ , , da 



r > r^ b < z < b + L 
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79 where a{r,z,t) is drawdown in the unsaturated zone, kQ{z) is relative permeabihty and 

80 Co{z) is moisture capacity (slope of the curve representing water saturation as a function 

81 of pressure head) functions with the functional dependence limitations on the respective 

82 constitutive models 

ko{z) = kieo), Co (z) = 0(60) (6) 

83 where Oq is the initial volumetric moisture content. Equation (|5| depends on the initial 

84 condition 

a(r,z,0) = 0, (7) 

85 the far-field boundary condition 

a{oo,z,t) = (8) 

86 the no-flow condition at the ground surface 

= r > (9) 

and the no-flow condition at the well casing 

b<z<b + L. (10) 
The interface conditions providing continuity across the water table are 



da 

dz 



z=h+L 



s — CT = r > r,,, z = b, 



'IV. 



ds da 



r > r,,, z = b. 



dz dz 

The aquitard drawdown si (r, z, t) is governed by 

I d f dsi\ d'^si „ 5si ^ ^ , ^ 

r dr \ dr I dz"^ dt 



(12) 



(13) 



90 Additionally, aquitard flow satisfies no-flow conditions at the bottom and center of the flow 

91 system 

= 0. (14) 



dsi\ dsi 
hm r— = — 
>o V dr J dz 



92 The interface condition across the aquifer-aquitard boundary are 

s — Si = r >ru, z = Q (15) 

93 and 

K,— = K,i-^ r>r^ z = Q. (16) 
oz oz 

94 Like Mishra and Neuman [5], we represent the aquifer moisture retention curve using an 

95 exponential function 



S^=^-M^ = e-.i^-^.) a^>o ^,>0 (17) 

by 

96 where 6r is residual volumetric water content, Sy = 6s — Or is drainable porosity or spe- 

97 cific yield and 5*6 is effective saturation. We also adopt the exponential relative hydraulic 

98 conductivity model [TO] . 



gafeW-i/'fc) ijj < 1p 



k 



Hip) ={ afc > V^fc > 0, 

1 il^> ipk 



99 with parameters and ipk that generally differ from Oc and ijja in (17). The parameters 

100 ak and Oc represent the exponent in the exponential models for hydraulic conductivity and 

101 effective saturation, respectively. The parameter ipk represents a pressure head above which 

102 relative hydraulic conductivity is effectively equal to unity, which is sometimes but not always 

103 equal to the air entry pressure head Va-In addition to rendering the resulting equations 

104 mathematically tractable, these exponential constitutive models are sufficiently flexible to 

105 provide acceptable fits to standard constitutive models such as those mentioned earlier. 

106 2.2. Point drawdown in saturated and unsaturated zones of aquifer and aquitard 

107 Following Mishra and Neuman [9], it is shown in Appnendix A that, drawdown in the 

108 saturated zone can be expressed as 

s = sc + su (19) 
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109 where sc is solution for flow to a partially penetrating well of finite radius in a confined aquifer 
no and su is a correction accounting for the underlying aquitard, water table and unsaturated 
111 zone effects. The Laplace transformed solution sc is given by Mishra and Neuman [9] as 



Sc [rD, zd,Pd) 



Q 



Ko(0o) + ^ Ci;„Ko(0„) cos [nn{l - Zd)] 



(20) 



n=l 



112 where 



sm{mxl£))—s\n(mTd£)) 



Vtin) 



113 2/Q(0), r^D = r^jr, rn = r/b, zd = z/b, = d/b, Id = l/b, po = pt, C^d = C^/iTrSsbrl), 

114 ts = ast/r"^, p is Laplace parameter (the transform of t), (f)n = \/pD/ts + rj^KDn'^Tr'^, and Kq 

115 and Ki are second-kind modified Bessel functions of orders zero and one. The Laplace trans- 

116 formed unsaturated zone drawdown a is given by Mishra and Neuman P] and is presented 

117 in Appendix D for sake of completeness. 

118 The Laplace transformed su derived in Appendix B is 



Su [J-d,zd,Pd) 



yK]i'^rD 



(21) 



119 where pi 



120 -R/<^/ii tanh (yUii?f,), 



A 

,,2 



P2 



, qib 



PD 



121 asi/as, Rt = bi/b, a^i = Kri/Ssi, and A = (|; + l) - l) e"^ - - 1 

122 The Laplace transformed aquitard drawdown derived in Appendix C is 



^ + 1 ) e^. 

qib I 



Sl[rD,ZD,PD) 



{Sc)zD=0 + Pl+ P2 

cosh(/ii6i/6) 



rhKn , 
X — ^yJo 



cosh [p,i{zD + Rb)] 



dy 



(22) 



123 where (sc)z^=o is the Laplace-Hankel transformed confined aquifer drawdown and is defined 

124 in Appendix D. The time domain equivalents sc, su, si and a of sc , su, si and a are 

125 obtained through numerical Laplace transform inversion using the algorithm of de Hoog 

126 et al. 



7 



127 2.3. Vertically Averaged Drawdown in Piezometer or Observation Well 

128 Drawdown in an observation well (Figure [T]) that is completed in the aquifer between 

129 elevations zoi = Zi/h and Z£i2 = Z2/h is found by averaging the point drawdown over screen 

130 interval, 

SzD2-ZDii'^'D,ts) = / s''{rD,ZD,ts) dZD (23) 

131 where s* can be either aquifer drawdown s, aquitard drawdown si, or a combination of the 

132 two, depending on the observation well screen interval. 

133 2.4. Delayed Piezometer or Observation Well Response 

134 When water level is measured in a piezometer or observation well having storage coeffi- 

135 cient C the water level observed in the borehole is delayed in time. Following Mishra and 

136 Neuman [9j, the measured (delayed) drawdown Sm can be expresses in terms of formation 

137 drawdown s via 

s^ = s(l-e-*/*^) (24) 

138 where is basic (characteristic) monitoring well time lag. The dimensionless equivalent of 



139 (|24j) is 

SmD = sd{1- e-*»/*^=) (25) 

140 where tss = ctstB/r"^, and r is the radial distance to the monitoring location. 



141 3. Model-predicted drawdown behavior 

142 We illustrate the impacts of an underlying aquitard on unconfined aquifer drawdown for 

143 the case where Kd = 1, Sgb/Sy = 10~^, akD = clcD = 10, ipao = "ipko, do = 0, C^jd = 10^, 

144 Id = 0.6 and = 0.02, where akD = cikb, a^D = CLcb, ipaD = 'ipa/b and ipkD = i^k/b. We also 

145 investigate the effects that wellbore storage capacity of the pumping well, the unconfined 

146 aquifer, and the unsaturated zone have on aquitard drawdown. 
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147 3.1. Dimensionless unconfined aquifer time-drawdown 

148 We start by considering drawdown at two locations in the unconfined aquifer saturated 

149 zone, one location closer to water table {z^, = 0.75) and the other closer to the aquitard- 

150 aquifer boundary (zd = 0.25). Figures [2^ and [2|d compare variations in dimensionless 

151 drawdown soiru, Zd,ts) = {4:TTKrb/Q)s{rD, ZD,ts) with dimensionless time at zd = 0.75 and 

152 Zd = 0.25 predicted by our proposed solution and the solutions of Mishra and Neuman 

153 [9], Neuman [27], and the modified solution of Malama et al. [21] (modified to include the 

154 partially penetrating pumping well effects, as done in Malama et al. [28] for a multi-aquifer 

155 system). The solutions of Neuman [2^ and Malama et al. [24J do not include wellbore storage 

156 effects, and therefore they overestimate drawdown at early time. Both of these solutions also 

157 ignore the unsaturated zone above the water table, considering the water table a material 

158 boundary [27]. Our proposed solution follows Mishra and Neuman [9] when leakage effects 

159 are minor, but our solution predicts less drawdown when leakage effects are significant. It 

160 is seen in Figure [2]d that solution of Mishra and Neuman [H] overestimates drawdown near 

161 the aquitard at intermediate time because it does not include aquitard leakage. Near the 

162 water table (Figure [2^) the effects of aquitard leakage are minimal and our proposed solution 

163 approaches Mishra and Neuman [9j at all times. 

164 Figures |3^ and show dimensionless time-drawdown variations at dimensionless radial 

165 distance ro = 0.5 and dimensionless unconfined aquifer saturated zone elevation zd = 0.25 

166 with different values of Rk^ = K^i/K^ when the radial aquitard hydraulic conductivity is 

167 small {Rxr = Kri/Kr = 10^^) and large {Rkr = 1-0). When the radial hydraulic conductivity 

168 in aquitard is negligible {Rkr = 10~^), aquitard flow is predominately vertical; larger values 

169 of vertical aquitard hydraulic conductivity cause decreases in intermediate time drawdown 

170 (Figure [3^). It is seen from Figure |3}d that when aquitard horizontal hydraulic conductivity is 

171 large {Rkr = 1) the amount drawdown is reduced from further increases in aquitard vertical 

172 hydraulic conductivity also extend to the later time. 

173 Figure |4] depicts the effect of Rxr on the dimensionless time-drawdown at dimensionless 

174 radial distance = 0.5 and dimensionless unconfined aquifer saturated zone elevation 
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175 Z£, = 0.25 when R^^ = 0.1. Radial flow in the aquitard results in less drawdown at late time 

176 than that predicted by Mishra and Neuman [9], who do not account for aquitard leakage. 

177 Figure |5] presents the effect of hydraulic conductivity of an isotropic aquitard on di- 
ns mensionless time-drawdown at dimensionless radial distance = 0.5 and dimensionless 

179 unconfined aquifer saturated zone elevation z^, = 0.25. When aquitard hydraulic conduc- 

180 tivity is at least two orders of magnitude smaller than the unconfined aquifer, the effects 

181 of leakage on the aquifer drawdown are negligible. This is in agreement with findings of 

182 Neuman and Witherspoon [16j for confined aquifers. They found errors < 5% attributable 

183 to the vertical aquitard fiow assumption, when the hydraulic conductivity contrast between 

184 the aquifer and aquitard was greater than a factor of 100. Figure |5] also presents a case 

185 with the aquitard hydraulic conductivity is larger than the aquifer. Because the proposed 

186 model accounts for general three-dimensional fiow in underlying zone, we can consider the 

187 case where the lower layer is more permeable than the aquifer {Rkr = 2). 

188 Figure |6] shows how the dimensionless unconfined aquifer time-drawdown is affected by 

189 aquitard thickness. When the aquitard thickness is less than the initial unconfined aquifer 

190 saturated thickness {Rb < 1) aquitard leakage only affects the time-drawdown curve at 

191 intermediate time. Figure [6] shows that further increases in aquitard thickness beyond eight 

192 times the initial unconfined aquifer saturated zone thickness have negligible effect on the 

193 time-drawdown curve. 

194 3.2. Dimensionless aquitard time- drawdown 

195 Figure[7]depicts dimensionless aquitard drawdown sniro, z^, tg) = {4:7iKrb/Q)si{rD, zd, ts) 

196 variations with dimensionless time at dimensionless radial distance td = 0.2 and dimension- 

197 less aquitard elevation zd = —0.25 for different values of C^d ■ As with solution of Mishra 

198 and Neuman |9J for non-leaky systems, aquitard drawdown is impacted by pumping-well 

199 wellbore storage capacity. Larger wellbore storage factors result in increased capacity of the 

200 wellbore to store water, resulting in a delay in the aquitard time-drawdown, as indicated in 

201 Figure 7. 

202 Figure [8] depicts the effect that changes in akD, the dimensionless relative hydraulic con- 
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203 ductivity exponent, have on dimensionless time-drawdown at dimensionless radial distance 

204 td = 0.2 and dimensionless aquitard elevation zd = —0.25. For larger values of SkD, the 

205 unsaturated zone hydraulic conductivity decreases more rapidly as pressure becomes more 

206 negative, relative to the threshold pressure ipk ■ A diminishing rate of water then drains from 

207 the vadose zone into the aquifer; this drainage contributes to reduced aquitard drawdown. 

208 For very large ako, unsaturated hydraulic conductivity quickly decreases once pressure head 

209 is below ipk , which leads to an much less permeable unsaturated zone. 

210 Figure |9] shows the effects that changes in a^D, the dimensionless effective saturation 

211 exponent, have on dimensionless time-drawdown at dimensionless radial distance td = 0.2 

212 and dimensionless aquitard elevation zd = —0.25 . When aco and ako are both large, 

213 pressure head and hydraulic conductivity in the vadose zone quickly reduce as pressure 

214 reaches the thresholds ipk and ipa- The vadose zone can no longer store water, and the water 

215 table essentially becomes a moving boundary, which leads to the limiting-case behavior of 

216 instantaneous drainage due to Neuman [27]. Consequently, for large values of exponents 

217 (Figure |9| red curve) the proposed solution reduces to that of Malama et al. [21], which 

218 relies on the assumption of instantaneous drainage of Neuman [27]. As decreases, the 

219 vadose zone has increased capacity to store water, which diminishes the water table response 

220 and aquifer drawdown increases, compared to that predicted by Malama et al. [21]. 

221 4. Conclusions 



222 Our work leads to the following major conclusions: 

223 1. A new analytical solution was developed for axially symmetric saturated- unsaturated 

224 three dimensional radial flow to a well with wellbore storage that partially penetrates 

225 the saturated zone of a compressible vertically anisotropic leaky-unconfined aquifer. 

226 The solution accounts for both radial and vertical flow in the unsaturated zone and 

227 the underlying aquitard. 

228 2. Because the solution considers three-dimensional radial flow in the aquitard, any prop- 

229 erties may be assigned to the aquitard, allowing the solution to also be used to simulate 
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230 leakage from underlying non-aquitard layers (e.g., an unscreened aquifer region with 

231 different hydraulic properties). 

232 3. Aquitard leakage can lead to significant departures from solutions that do not account 

233 for leakage, e.g., Mishra and Neuman [U]. However, the effect of leakage on uncon- 

234 fined aquifer drawdown diminishes at points farther away from the aquifer-aquitard 

235 boundary. 

236 4. Unsaturated zone effects are often more important than leakage effects when the ob- 

237 servation location is close to the water table. 

238 5. For large diameter pumping wells, at early time water is withdrawn entirely from the 

239 wellbore storage. Solution that do not account for wellbore storage predict a much 

240 larger early rise in drawdown. 

241 6. Aquitard drawdown is also affected by the pumping- well wellbore storage capacity. As 

242 in the unconfined aquifer, larger wellbore storage capacity leads to larger impacts on 

243 the observed aquitard drawdown. 

244 7. The unsaturated zone properties not only affect the unconfined aquifer time-drawdown 

245 behavior but they also impact the observed aquitard response. 
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Pumping well 




Figure 1: Schematic representation of leaky unconfined aquifer-aquitard system geometry with finite radius 
pumping well 



255 Appendix A. Decomposition of saturated zone solution 

256 In a manner analogous to Mishra and Neuman [5] we decompose s into two parts 



s = Sc + Su 



(A.l) 



257 where sc is solution for a partially penetrating well in a confined aquifer, satisfying 

1 d ( dsc\ d'^sc I dsc n ^ / a 

r— — I + Kb ^ ^ = — r > r,„ < z < o 



r dr \ dr I dz^ as dt 

Sc{r, z,0) = r > 



(A.2) 
(A.3) 



Sc {oo, z,t) = 



ds 



c 



dz 



r > r„ 



z={0,b) 



(A.4) 
(A.5) 



dsc 
dr 



0< z <b-l b-d< z<b 



(A.6) 
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Table 1: Fundamental Properties Table 



a 


Hankel transform parameter 






exponent in moisture retention curve or sorptive number 


L-i 




exponent in Gardner relative hydraulic conductivity model 


L-i 


b 


saturated thickness of unconfined aquifer before pumping begins 


L 


h 


thickness of aquitard 


L 




wellbore storage coefficient 




d 


distance from top of screened interval to top of aquifer 


L 


h 


hydraulic head (sum of pressure and elevation heads) 


L 




aquifer radial hydraulic conductivity 


LT-^ 


Krl 


aquitard radial hydraulic conductivity 


LT-^ 




aquifer vertical hydraulic conductivity 


LT-^ 




aquitard vertical hydraulic conductivity 


LT-^ 


I 


distance from bottom of screened interval to top of aquifer 


L 


L 


thickness of vadose zone before pumping begins 


L 


n 


finite cosine transform parameter 


— 


P 


Laplace transform parameter 




Q 


volumetric pumping rate 




r 


radial distance from the center of pumping well 


L 




diameter of pumping well 


L 


S 


drawdown in aquifer; change in hydraulic head since pumping began 


L 


Si 


drawdown in aquitard; change in hydraulic head since pumping began 


L 


Se 


effective saturation 


— 


Ss 


aquifer specific storage 




Ssl 


aquitard specific storage 




Sy 


aquifer drainable porosity or specific yield 


— 


t 


time since pumping began 


T 


z 


vertical distance from the bottom of the aquifer, positive up 


L 




elevation to top {i = 1) and bottom {i = 2) of monitoring interval 


L 




initial volumetric water content 


- 


Of 


residual volumetric water content 


— 


Os 


saturated volumetric water content 




a 


drawdown in unsaturated zone; change in hydraulic head since pumping began 


L 




pressure head (less than zero when unsaturated) 


L 




air-entry pressure 


L 


'A: 


pr(>ssur(' Ibr saturated hydraulic conduct i\ it,y| ^ 


L 



Table 2: Derived quantities table 



Kd 


KjKr 


Anisotropy ratio 


ro 


r/b 


dimensionless radial coordinate 


ZD 


z/b 


dimensionless vertical coordinate 


do 


d/b 


dimensionless distance to top of screen interval 


Id 


l/b 


dimensionless distance to bottom of screen interval 


Pd 


pt 


dimensionless Laplace parameter 




rw/r 


dimensionless well radius 


Rro 


Kdi/Kd 


ratio of aquitard and aquifer anistropies 




Krl/Kr 


ratio of aquitard and aquifer horizontal hydraulic conductivities 




K^i/K, 


ratio of aquitard and aquifer vertical hydraulic conductivities 


Ras 




ratio of aquitard and aquifer saturated hydraulic diffusivities 


Rb 


bi/b 


ratio of aquitard and aquifer thicknesses 




Kr/Ss 


aquifer hydraulic diffusivity 


(Xsl 


Krl/Ssl 


aquitard hydraulic diffusivity 


^D. 


Zi/b 


dimensionless elevation to top (i = 1) and bottom {i = 2) of monitoring interval 


O-kD 


ttkb 


dimensionless Gardner hydraulic conductivity model exponent 


dcD 


aj) 


dimensionless moisture retention model exponent 


-tpaD 




dimensionless air-entry pressure 


tpkD 


tpk/b 


dimensionless pressure for saturated hydraulic conductivity 


CwD 


Cw/ijrSMl) 


dimensionless wellbore storage coefficient 


ts 




dimensionless time 
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10^ 




Dimensionless Time [ tg ] 




Dimensionless Time [ tg ] 



Figure 2: Dimensionless leaky-unconfined aquifer drawdown versus dimensionless time at rjj — 0.5 when 
Kd - 1, Ssb/Sy = 10-^ akD = a,D = 10, V^d = V^feD, do ^0,Id = 0.6, Cy,D = 10^ Rk,. = Rk^ = lO^^, 
Rs^ — lO^'^, i?h ^ oo and (a) zd — 0.75 (b) zo — 0.25. Solutions of Mishra and Neunian [9 , modified 
Malama et al. [24=, and Neuman [27] are also shown. 




Mishra and Neuman [2011] ■ 

Rkz =10"' 
Rkz =10 ■ 
Rk; = 0.50 
Rk, = 1.00 



10' 10 lO-' 10 

Dimensionless Time [ tg ] 



10= 









Misiira and Neuman [201 1] 
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Figure 3: Dimensionless leaky-unconfined aquifer drawdown versus dimensionless time at = 0.5 and 
ZD = 0.25 for Kd = 1, Ssb/Sy = 10~^, ako = Ucd = 10, ^paD = ipko, dD ^ 0, Id ^ 0.6, C^d = 10^, 
Rs^ = 10^^, i?f, — > oo when Rk, varies and (a) Rkr — 10^^ (b) Rk^ = 1. Solution of Mishra and Neuman 
[2] is also shown. 
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Figure 4: Dimensionless leaky-unconfined aquifer drawdown versus dimensionless time at ro = 0.5 and 
ZD = 0.25 for Kd = 1, Ssb/Sy = 10~^, ako = OcD = 10, tpaO = V'feD, do = 0, /d = 0.6, C'y,D = 10^, 
i?5^ = 10^^, Rk^ =0.1 and i?^ — > oo when Rx^ varies. Solution of Mishra and Neuman [9" is also shown. 




Figure 5: Dimensionless leaky-unconfined aquifer drawdown versus dimensionless time at = 0.5 and 
ZD = 0.25 for Kd = 1, Ssb/Sy = 10"^, ako = clcD = 10, i/'aD = V^feD, c?_d = 0, /d = 0.6, C„_d = 10^, 
i?b — > cx) when = -Rx^ varies and i?5^ = 1. Solution of Mishra and Neuman [9_ is also shown. 
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Figure 6: Dimensionless leaky-unconfined aquifer drawdown versus dimensionless time at ro ~ 0.5 and 
ZD = 0.25 for Kd = 1, Ssb/Sy = 10~^, ako = Ocd = 10, ipaD = V'feD, do = 0, /d = 0.6, C'y,D = 10^, 
Rs^ = 10^, Rk^ — Rk^ = 10^^ when Ri, = bi/b varies. Solution of Mishra and Neuman [91 is also shown. 




Figure 7: Dimensionless aquitard drawdown versus dimensionless time at ro — 0.2 and zjj = —0.25 for 
Kd = 1, Ssb/Sy = 10-3, akD = a,D = 10, xjjao - ^ko, do ^ 0, h = 0.6, i?s. = 10^ Rk^ = Rk,. = 10"', 
i?b — > cx) when C^d, the dimensionless wellbore storage varies. 
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Dimensionless Time [ ] 



Figure 8: Dimensionless aquitard drawdown versus dimensionless time at ro — 0.2 and zd = —0.25 for 
Kd = 1, Ssb/Sy = 10'^ ^PaD = i'kD, do ^0,Id = 0.6, C^d = 10'\ = 10^, Rk^ - i?AV = 10-^ 
— > oo when, = 1 and ako varies. 




Dimensionless Time [ tg ] 

Figure 9: Dimensionless aquitard drawdown versus dimensionless time at ro — 0.2 and zjj — —0.25 for 
Kd - 1, Ssb/Sy = 10-3, VaD = i'kD, do ^ 0, Id = 0.6, C^d = 10^, i?s, = 10^, Rk^ = Rk,. = 
— > oo when, a^jj = 10"^ and Ucd varies. 
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262 



2T^Kr {I - d) 



dsc 
dr 



ds 



c 



w ■ 



dt 



-Q 



b-l< z <b-d 



(A.7) 



263 and su is a solution that takes into account aquitard and saturated-unsaturated unconfined 

264 conditions, but has no pumping source term, satisfying 



1 (9 / ds 



'1/ 



r dr \ dr 



K 



d'^su 1 ds 



u 



D 



dz^ as dt 

r>0 Q< z<h 



s^(r,z,0) = r>0 



(A.9) 



Su (oo, z,t) — Q 

dsu _ K^i dsi 
dz Kz dz 

dsu 



=0 r>0 z^O 



dz 



< ^ < 6 



r=0 



269 subject to interface conditions at water table. 



sc + Su — cr = r > rw z = b 

dsc , dsu da 
dz dz dz 



(A. 10) 
(A.ll) 

(A.12) 



(A.13) 
(A.14) 



271 where the first term is zero by definition of sc- 
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272 Appendix B. Laplace-space solution for saturated zone 



273 Equations (A.8)-(A.14) are solved by sequential application of the Hankel transform 

(B.i: 



and Laplace transform 



/(a) = J rJo(ar)/(r) dr 





(B.2) 



fip) = J fit)e-^' dt 



275 with Hankel parameter a and Laplace parameter p, Jq being zero-order Bessel function of 

276 the first kind. 

277 The Laplace-Hankel transform of confined aquifer solution [S] is 



Sc[a,ZD,PD) 



Co \ — Ji(ar^)Ko(r^ro) 

ror^Jo(ar)Ki(rro) - ar^Ji(ar)Ko(rro) ) 

J 



a2 + r2 



CXJ 

Cn \ — Ji(ar^)Ko(r^ro) 
^ — ^ L a 



n=l 



r„r^Jo(ar)Ki(rr„) - ar^Ji(ar)Ko(rr„) 



a2 + r2 



X cos [mx{l — zd)] 



278 where tq = ^/pS^JK^ and Tn = \JpSs/Kr + Kj^v?"n'^ jh"^ 

279 The Laplace transform of (A.8)-(A.14) is 

d'^su p _ 



r dr \ dr 



su < z <b 



(B.3) 



(B.4) 



Su {oo, z,p) = 



(B.5) 



ds 



u 



dz 



z=0 



Kzi dsi 
dz 



1 =0 < z < 6 



dr 



r=0 



(B.6) 
(B.7) 
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sc + su — <y = z = b 
dsc dsu da 



dz dz dz 

285 where the first term is zero by definition of sc- 



z = b 



286 Taking the Hankel transform of (|B.4|)-(|B.9|) yields 

-Bu < z <b 



2= , d^^U P = 



- a su + Kd 



dz"" 



a. 



sh + Su — ci = Q z = b 
dsH dsu da 



dz ~'~ dz dz 

d^c _^ dsij_ _ K^i dsi 
dz dz K, dz 



z = b 



z = 



The general solution of (B.IO) subject to (B.ll) is 



Su = Pie' + P2e 



-rjz 



291 where pi and p2 are coefficients to be determined from boundary conditions. 



292 Considering that dsu/dz = at 2; = and z = bhj virtue of (A.5) and that 



dz 



z=b 



293 which, together with g, are derived in (D15) of Mishra and Neuman |9j and 

dMu 



dz 



qi[sc + Su, 



2=0 



z=0 



294 which, together with qi , are derived in (C.7) we obtain from (B.ll )-(B.13) 

Qiiv + q)e~''^sc{z = 0) - q{r] + qi)s{z = b) 



Pi 



P2 



A 



Qiiv - g)e ''^ scjz = 0) - q{ri - qi)s{z = b) 
A 
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296 where A = {t] — qi){rj + q)e ^'^ — {j] — q){rj + qi)e'^^. 

297 The inverse Hankel transform of ( |B.14 ) is 

(pie'''' + P2e~'^'') aJo(ar) da. 



su 



(B.19) 



1/2 I I I 

298 Defining a new variable y = ar/K^ ro transforms ( |B.19 ) into the resuh presented in (21). 

299 It is noted that when gi = the aquiard is replaced by an impermeable boundary, and 



300 Pi = P2 



2sc{z=h) 
cosh(?76) — - sinh(»jfe) 



. These simplifications reduce (B.19) to equation (3) of Mishra 



and Neuman [9]. 

Appendix C. Aquitard Solution 

Laplace-Hankel transform of governing flow equations for aquitard are 

P 

By virtue of no flow boundary at the bottom of the system, ^|^__^ = 0, the 



-a si + Km—— = p-—si < z < -6i 



305 solution to (C.l) is 



where ril 



si = pi cosh [?7i {z + hi)] 
-w— + J*^^ ■ The boundary condition 

tl{z = 0) = 1(2; = 0) = (Ic + 1(7)^=0 



307 gives 



cosh(?76i) 



cosh \r]i{z + &i)] . 



(C.l) 
general 

(C.2) 

(C.3) 
(C.4) 



Using (B.13) transforms (C.4) into the solution presented in (22). 



The derivative of (C.4) is 

dsi 
dz 



?7itanh(r7i&i)(sc + su) 



z = 0. 



The flux boundary condition at the aquifer-aquitard interface 



dsi 

dz 



ds 



2=0 



2=0 



Combined with (B.9) gives 



ds 
dz 



Kzi dz 



(C.5) 



(C.6) 



(C.7) 



2=0 



312 where qi = ^771 tanh (7716 
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313 Appendix D. Lapalce Transformed Unsaturated Zone Drawdown 



314 Lapalce transformed drawwdown cr in tlie unsaturated zone are given by Mishra and 

315 Neuman [9] as 



o-{rD,ZD,PD) 



/c 



°° ^akn(zn-l)/2 ^n[i(p{zD-l)]+xyn[i<f'{zD-'i-)] 



Jn[i</.(0)]+xY„[i</.(0)] 



X {sc + Su) 



2_D=1 



dy 



for ^ CLkD 



(D.i: 



/c 



oo e''lD(^-D~l)+xe''2D(^-D-l) 
1+X 



316 where rr, = r/b, zn = z/b, /i^ = v"^ + ^ 

317 akD = akb, a^D = cicb, 4>{zd) = 



ast/r'^, as = Kr/Ss, go = qb 



VP y- 



tsKor' 



318 Sd = Sy/S, IpkD = ^k/b, IpaD = -ipa/b, 5lD,2D 

319 and Pd = pt are dimensionless quantities, p being the Laplace transform parameter; 

„2 



ic(z = b) = a 



r 



KdtI 



Ji {yDTwD) Ko (r^„D0o) + r(o) 



n=l 



y2 



Ji {yDTwD) Ko (r^D0o) + r(ra) 



(D.2) 



where r(n) 



'rwD<t>nia(yD)'^l{4>n)-yDr'wDil(yD)^0{4'n) 



Vd 



1 /2 

r£), Jo and Ji being Bessel 



functions of first kind and, respectively, orders zero and one; and 



(afcO+nAD)Jn[^(/>(Lp)]-2^A/Boe^-P^oJ„+l[t(/.(Lp)] 
{akD+n-XoWn [i4>{L]j)]- 



x= < 



akD 7^ o-cD, Ld ^ 00 
akD = acD 

akD = acD, Ld — > 00 



(D.3) 



Qd 



( akD I nXD \ _ • /Q— Jn+i[''/'(0)]+xYn+l[»'j!'(0)] / 
{2^2) ''V^D J„[j0(O)]+xY„[i<^(O)] akDTacD 



1+X 



(D.4) 



akD — acD — i^D 

322 where Ld = L/b, Jn and Y„ being first and second kind Bessel functions of order n. 
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